clear; clc; 

%% Fig 1
% panel (A)
para00 = [1 1];
N_mkt = 4000;
N_vec = (2:4);
phi_vec = [0.4,0.5,0.6,0.7];
N_monte = 300;

figure; set(gcf, 'Position',  [100, 100, 1100, 900], 'PaperUnits', 'inches','PaperSize', [15 10]);
tiledlayout(length(N_vec), length(phi_vec), "TileSpacing","compact");
for j2 = 1:length(N_vec)
    for j1 = 1:length(phi_vec)
        nexttile;
        
        tmp = sprintf('MonteCarloEstResult_N%1.0f_Nmonte%1.0f_Nmkt%1.0f_phi%1.2f.mat', N_vec(j2), N_monte, N_mkt, phi_vec(j1));
        load(tmp);

        CoverageProb_fy = nan(length(Xgrid), 1);
        CoverageProb_ct = nan(length(Xgrid), 1);
        CoverageProb_at = nan(length(Xgrid), 1);
        for i = 1:length(Xgrid)
            tmp = any(covered_mesh_fy_as10(:, i, :), 1);
            CoverageProb_fy(i) = squeeze(mean(tmp, 3));
        
            tmp = any(covered_mesh_ct_as10(:, i, :), 1);
            CoverageProb_ct(i) = squeeze(mean(tmp, 3));
        
            tmp = any(covered_mesh_at_as10(:, i, :), 1);
            CoverageProb_at(i) = squeeze(mean(tmp, 3));
        end
        plot(Xgrid(:) - para00(1), CoverageProb_fy(:), 'b-', 'LineWidth',2); hold on;
        plot(Xgrid(:) - para00(1), CoverageProb_ct(:), 'k--', 'LineWidth',2); hold on;
        plot(Xgrid(:) - para00(1), CoverageProb_at(:), 'r:', 'LineWidth',2); hold off;
        legend('FY','CT','AT', 'FontSize',8);
        xlabel('$C - C_{true}$','Interpreter','latex');
        ylabel('Coverage Probability');
        title(sprintf('$N = %1.0f, \\phi = %1.1f$', N_vec(j2), phi_vec(j1)), 'Interpreter','latex', 'FontSize',16);
        pbaspect([1 0.8 1]);
    end
end
print('CoverageProb_FcMean.pdf', '-dpdf','-bestfit');

% Panel (B)
figure; set(gcf, 'Position',  [100, 100, 1100, 900], 'PaperUnits', 'inches','PaperSize', [15 10]);
tiledlayout(length(N_vec), length(phi_vec), "TileSpacing","compact");
for j2 = 1:length(N_vec)
    for j1 = 1:length(phi_vec)
        nexttile;
    
        tmp = sprintf('MonteCarloEstResult_N%1.0f_Nmonte%1.0f_Nmkt%1.0f_phi%1.2f.mat', N_vec(j2), N_monte, N_mkt, phi_vec(j1));
        load(tmp);
        
        CoverageProb_fy = nan(length(Ygrid), 1);
        CoverageProb_ct = nan(length(Ygrid), 1);
        CoverageProb_at = nan(length(Ygrid), 1);
        for i = 1:length(Ygrid)
            tmp = any(covered_mesh_fy_as10(i, :, :), 2);
            CoverageProb_fy(i) = squeeze(mean(tmp, 3));
        
            tmp = any(covered_mesh_ct_as10(i, :, :), 2);
            CoverageProb_ct(i) = squeeze(mean(tmp, 3));
        
            tmp = any(covered_mesh_at_as10(i, :, :), 2);
            CoverageProb_at(i) = squeeze(mean(tmp, 3));
        end
        plot(exp(Ygrid(:)) - para00(1), CoverageProb_fy(:), 'b-', 'LineWidth',2); hold on;
        plot(exp(Ygrid(:)) - para00(1), CoverageProb_ct(:), 'k--', 'LineWidth',2); hold on;
        plot(exp(Ygrid(:)) - para00(1), CoverageProb_at(:), 'r:', 'LineWidth',2); hold off;
        xlabel('$\sigma - \sigma_{true}$','Interpreter','latex');
        legend('FY','CT','AT', 'FontSize',6);        
        ylabel('Coverage Probability');
        title(sprintf('$N = %1.0f, \\phi = %1.1f$', N_vec(j2), phi_vec(j1)), 'Interpreter','latex', 'FontSize',16);
        pbaspect([1 0.8 1])
    end
end
print('CoverageProb_FcStd.pdf', '-dpdf','-bestfit');

%% Fig 2
para00 = [1 1 1];
N_mkt = 4000;
phi = 0.5;
N = 2;
N_monte = 100;

% Left panel
figure; 
tmp = sprintf('MonteCarloEstResult_N%1.0f_Nmonte%1.0f_Nmkt%1.0f_phi%1.2f.mat', N, N_monte, N_mkt, phi);
load(tmp);
CoverageProb_fy = nan(length(Xgrid), 1);
CoverageProb_ct = nan(length(Xgrid), 1);
CoverageProb_at = nan(length(Xgrid), 1);
for i = 1:length(Xgrid)
    tmp = any(any(covered_mesh_fy_as10(:, i, :, :), 1), 3);
    CoverageProb_fy(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_ct_as10(:, i, :, :), 1), 3);
    CoverageProb_ct(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_at_as10(:, i, :, :), 1), 3);
    CoverageProb_at(i) = squeeze(mean(tmp, 4));
end
plot(Xgrid(:) - para00(1), CoverageProb_fy(:), 'b-', 'LineWidth',2); hold on;
plot(Xgrid(:) - para00(1), CoverageProb_ct(:), 'k--', 'LineWidth',2); hold on;
plot(Xgrid(:) - para00(1), CoverageProb_at(:), 'r:', 'LineWidth',2); hold off;    
legend('FY','CT','AT', 'FontSize',6);
xlabel('$C - C_{true}$','Interpreter','latex');
ylabel('Coverage Probability');
pbaspect([1 1 1])
print(['CoverageProb_FcMean_corrN' num2str(N) '.pdf'], '-dpdf');

% middle panel
figure; 
CoverageProb_fy = nan(length(Ygrid), 1);
CoverageProb_ct = nan(length(Ygrid), 1);
CoverageProb_at = nan(length(Ygrid), 1);
for i = 1:length(Ygrid)
    tmp = any(any(covered_mesh_fy_as10(i, :, :, :), 2),3);
    CoverageProb_fy(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_ct_as10(i, :, :, :), 2),3);
    CoverageProb_ct(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_at_as10(i, :, :, :), 2),3);
    CoverageProb_at(i) = squeeze(mean(tmp, 4));
end
plot(exp(Ygrid(:)) - para00(2), CoverageProb_fy(:), 'b-', 'LineWidth',2); hold on;
plot(exp(Ygrid(:)) - para00(2), CoverageProb_ct(:), 'k--', 'LineWidth',2); hold on;
plot(exp(Ygrid(:)) - para00(2), CoverageProb_at(:), 'r:', 'LineWidth',2); hold off;
xlabel('$\sigma - \sigma_{true}$','Interpreter','latex');
legend('FY','CT','AT', 'FontSize',6);
ylabel('Coverage Probability');
pbaspect([1 1 1]);
print(['CoverageProb_FcStd_corrN' num2str(N) '.pdf'], '-dpdf');

% right panel
figure; 
CoverageProb_fy = nan(length(Zgrid), 1);
CoverageProb_ct = nan(length(Zgrid), 1);
CoverageProb_at = nan(length(Zgrid), 1);
for i = 1:length(Zgrid)
    tmp = any(any(covered_mesh_fy_as10(:, :, i, :),1), 2);
    CoverageProb_fy(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_ct_as10(:, :, i, :), 1), 2);
    CoverageProb_ct(i) = squeeze(mean(tmp, 4));

    tmp = any(any(covered_mesh_at_as10(:, :, i, :), 1), 2);
    CoverageProb_at(i) = squeeze(mean(tmp, 4));
end
plot(exp(Zgrid(:)) - para00(3), CoverageProb_fy(:), 'b-', 'LineWidth',2); hold on;
plot(exp(Zgrid(:)) - para00(3), CoverageProb_ct(:), 'k--', 'LineWidth',2); hold on;
plot(exp(Zgrid(:)) - para00(3), CoverageProb_at(:), 'r:', 'LineWidth',2); hold off;
xlabel('$\sigma^{mkt} - \sigma^{mkt}_{true}$','Interpreter','latex');
legend('FY','CT','AT', 'FontSize',6);
ylabel('Coverage Probability');
pbaspect([1 1 1])
print(['CoverageProb_FcMktStd_corrN' num2str(N) '.pdf'], '-dpdf');

